clear all
close all
%% RPM traindata
% Generate training dataset using rock physics model
% define ROCK PHYSICS INPUTS
ntrain=801; % number of points
sghsim = rand(ntrain,1);
phisim = rand(ntrain,1).*0.6;
claysim = ones(ntrain,1).*0.6;
phic = 0.4; % critical porosity
n = 4; % coordination number
peff =  rand(ntrain,1)*0.005; % 0-5 MPa
% ROCK PHYSICS FORWARD MODEL
clear rhosim vpsim vssim
[rhosim,vpsim,vssim]=Dvorkin_FT_loadbearing(phisim,phic,claysim,sghsim,peff,n);
modsim = [phisim sghsim claysim];
datasim = [vpsim, vssim, rhosim];
traindata = [modsim datasim];  % this is the traning dataset (created using the RPM)
traindata(isnan(traindata(:)) == 1) = 0;
ncomp = 4; % for Gaussian Mixture Model, Grana (2016)
%% load realizations
ProjectPath = 'G:\backup_win_d\ava_inversion\apb17\';
InputPath   = 'input\';
OutputPath  = 'output_RHO_HD_LFM\';

AvgVp = zeros(1,1613,801);
AvgVs = zeros(1,1613,801);
AvgDen = zeros(1,1613,801);

for i=1:32
    vp = sgems_read([ProjectPath OutputPath 'vp_DSS_6_',num2str(i),'_1.sgems']);
    vs = sgems_read([ProjectPath OutputPath 'vs_coDSS_6_',num2str(i),'_1.sgems']);
%     den = sgems_read([ProjectPath OutputPath 'd_coDSS_6_',num2str(i),'_1.sgems']);
    den = sgems_read([ProjectPath OutputPath 'best_D_2.sgems']); % fixed RHO
  
    iter6(i).vp = reshape(vp.data,1,1613,801);
    iter6(i).vs = reshape(vs.data,1,1613,801);
    iter6(i).den = reshape(den.data,1,1613,801);
    % Trace by trace inversion
    clear *section
    for j=1:size(iter6(i).den,2)
        clear obsdata phiest phiml sghest sghml
        obsdata = [iter6(i).vp(1,j,:) iter6(i).vs(1,j,:) iter6(i).den(1,j,:)];
        obsdata = squeeze(obsdata)';
        i
        j
        [phiest, ~, sghest, ~,clayss] = HydrateInversion(obsdata,traindata, ncomp);
        phisection(1,j,:)    = phiest;
        sghsection(1,j,:)    = sghest;
    end
    
    iter6(i).phi = phisection;
    iter6(i).sgh = sghsection;
    
    AvgVp = AvgVp + iter6(i).vp;
    AvgVs = AvgVs + iter6(i).vs;
    AvgDen = AvgDen + iter6(i).den;
end

AvgVp = AvgVp./32;
AvgVs = AvgVs./32;
AvgDen = AvgDen./32;

% save('G:\backup_win_d\ava_inversion\OK_volumetrics_21_17\UQ_apb21_realiz.mat','iter6','-v7.3')
save('G:\backup_win_d\ava_inversion\OK_volumetrics_21_17\UQ_apb17_realiz_fixedRHO.mat','iter6','-v7.3')
